function [params] = Newton_Raphson(myfunc, params0, tolX, iterMax, varargin)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% This function applies the Newton-Raphson method to optimize a function,
	% whose gradient and hessian need to be calculated.
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Inputs:
	% myfunc:				function:
	%							takes		NumParams x 1 as input
	%							returns:	objective:	scalar
	%										gradient:	NumParams x 1
	%										hessian:	NumParams x NumParams
	% params0:				NumParams x 1
	% tolX:					scalar
	% iterMax:				integer
	% varargin{1}=stepSize:	scalar (between 0 and 1)
	% varargin{2}=verbose:	boolean
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Outputs:
	% params:				NumParams x 1
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	% Read optional arguments
	if length(varargin) >= 1; stepSize = varargin{1}; else; stepSize = 1;    end
	if length(varargin) >= 2; verbose  = varargin{2}; else; verbose = false; end
	
	% Do initializations
	params = params0;
	iter = 0;
	criterionX = Inf;
	
	while iter < iterMax && criterionX > tolX
		[obj, grad, hess] = myfunc(params);
		params_new = params - stepSize * hess\grad;
		iter = iter + 1;
		criterionX = max(abs(params_new-params));
		params = params_new;
		if verbose
			disp(sprintf('Newton-Raphson iteration %d: obj=%g, criterion=%g', iter, obj, criterionX));
		end
	end
	disp(sprintf('Newton Raphson algorithm has finished afer %d iterations.', iter));
	disp(sprintf('Value of objective: %g', obj));
	disp(sprintf('Value of criterion: %g', criterionX));
end
